Modified pulse sequence to estimate properties

ABSTRACT

Methods and related systems are described for estimating fluid or rock properties from NMR measurements. A modified pulse sequence is provided that can directly provide moments of relaxation-time or diffusion distributions. This pulse sequence can be adapted to the desired moment of relaxation-time or diffusion coefficient. The data from this pulse sequence provides direct estimates of fluid properties such as average chain length and viscosity of a hydrocarbon. In comparison to the uniformly-spaced pulse sequence, these pulse sequences are faster and have a lower error bar in computing the fluid properties.

PRIORITY

This application claims priority as a continuation application of U.S.patent application Ser. No. 13/543,393, filed Jul. 6, 2012, which is acontinuation of Ser. No. 12/717,512, filed Mar. 4, 2010, now granted asU.S. Pat. No. 8,587,302.

BACKGROUND

1. Field

This patent specification relates to estimating properties fromexponentially decaying data. More particularly, this patentspecification relates to methods and systems for estimating materialproperties from a non-uniform sequence of excitations, such as with NMRmeasurements.

2. Background

Traditionally, fluid properties using Nuclear Magnetic Resonance (NMR)are obtained by study of longitudinal or transverse relaxation ordiffusion measurements. Longitudinal (T₁) relaxation data are acquiredusing an inversion-recovery pulse sequence where the recovery times areoften chosen arbitrarily. Transverse (T₂) relaxation data are acquiredusing the CPMG pulse sequence with pulses that are typically uniformlyspaced in the time-domain. There is no clear mathematical methodologyfor choosing the acquisition times at which the diffusion measurementsare obtained. The relaxation or diffusion data thus obtained from thesepulse sequences are used to infer the relaxation time or diffusiondistributions. In turn, these distributions are used to inferpetrophysical or hydrocarbon properties.

Traditional well-logging using NMR is employed to infer petro-physicaland fluid properties of the formation or in-situ hydrocarbon. Themulti-exponential time-decay of NMR magnetization is characterized byrelaxation time constants T₁, T₂ or diffusion D and correspondingamplitudes ƒ_(T) ₁ (T₁), ƒ_(T) ₂ (T₂) or ƒ_(D)(D). Since the relaxationtimes and/or diffusion constants are expected to be continuous andtypically span several decades, these amplitudes are often referred toas relaxation-time or diffusion distributions. These distributionsprovide a wealth of information about both the rock as well as fluidproperties. See e.g. D. Allen, S. Crary, and R. Freedman, How to useborehole Nuclear Magnetic Resonance, Oil eld Review, pages 34-57, 1997.

An NMR measurement consists of a sequence of transverse magnetic pulsestransmitted by an antenna. Typically, the CPMG pulse sequence is used tomeasure T₂ relaxation. It consists of a pulse that tips hydrogen protons90° and is followed by several thousand pulses that refocus the protonsby tipping them 180°. The data are acquired by an antenna between theevenly spaced 180° pulses and are thus uniformly spaced in time. On theother hand, the T₁ relaxation is studied using the inversion-recoverypulse sequence. There is no specific rule for choosing the recoverytimes. See, Y. Q. Song, L. Venkataramanan, and L. Burcaw, Determiningthe resolution of Laplace inversion spectrum, Journal of ChemicalPhysics, page 104104, 2005. Similarly, there is no clear mathematicalframework to optimally choose the acquisition times for diffusionmeasurements.

The data acquired from all of the above pulse sequences can berepresented by a Fredholm integral (see, M. D. Hurlimann and L.Venkataramanan, Quantitative measurement of two dimensional distributionfunctions of di usion and relaxation in grossly homogeneous elds,Journal Magnetic Resonance, 157:31-42, July 2002 and L. Venkataramanan,Y. Q. Song, and M. D. Hurlimann, Solving Fredholm integrals of the rstkind with tensor product structurein 2 and 2.5 dimensions, IEEETransactions on Signal Processing, 50:1017-1026, 2002, hereinafter“Venkataramanan et al. 2002”),

$\begin{matrix}{{M(t)} = {\int_{0}^{\infty}{e^{{- t}/x}{f_{x}(x)}{dx}}}} & (1)\end{matrix}$

where the variable x typically refers to T₁, T₂, or 1/D. See,Venkataramanan et al. 2002. Thus the measured data in eqn. (1) isrelated to a Laplace transform of the unknown underlying functionƒ_(x)(X).

Scaling laws established on mixtures of alkanes provide a relationbetween fluid composition and relaxation time and/or diffusioncoefficients of the components of the mixture (see D. E. Freed, Scalinglaws for di usion coe cients in mixtures of alkanes, Physical ReviewLetters, 94:067602, 2005, and D. E. Freed, Dependence on chain length ofNMR relaxation times in mixtures of alkanes, Journal of ChemicalPhysics, 126:174502, 2007). A typical crude oil has a number ofcomponents, such as methane, ethane etc. Let N_(i) refer to the numberof carbon atoms or chain length of the i-th component. Let Ndenote themolar average chain length of the crude oil, also defined as theharmonic mean of the chain lengths,

$\begin{matrix}{{\frac{1}{\overset{\_}{N}} \equiv {\sum\limits_{i = 1}^{\infty}\frac{f_{N}\left( N_{i} \right)}{N_{i}}}} = \left\langle N^{- 1} \right\rangle} & (2)\end{matrix}$where ƒ_(N)(N_(i)) denotes the mass fraction of component i in thehydrocarbon. According to the scaling law, the bulk relaxation timeT_(2,i) of component i in a crude oil follows a power law and dependsinversely on its chain length or number of carbon atoms, N_(i). Thus,the larger the molecule, the smaller the relaxation time T_(2,i). Therelaxation time also depends inversely on the average chain length N ofthe hydrocarbon, according toT _(2,i) =BN _(i) ^(−κ) N ^(−γ)  (3)

In eqn. (3), parameters B, γ and κ are constants at a given pressure andtemperature. Similarly, if D_(i) refers to the diffusion coefficient ofcomponent i,D _(i) =AN _(i) ^(−ν) N ^(−β)  (4)where A, ν and β and are constants at a given pressure and temperature.The scaling theory can also be used to derive an expression for fluidviscosity,η∝

D ⁻¹

η∝N ^(β)

N ^(Γ)

.  (5)

Traditional analysis of the measured NMR data to obtain fluid propertiesgoes as follows. The inverse Laplace transform of eqn. (1) is performedto estimate the probability density function η_(x)(x) from the measureddata. Next, using eqns. (3)-(5), the average chain length N and thechain length distribution ƒ_(N)(N) are computed. Viscosity η of thehydrocarbon is then estimated using eqn. (5). This traditional analysishas some disadvantages. First, the inverse Laplace transform of themeasured NMR data in eqn. (1) is non-linear due to the non-negativityconstraint on ƒ_(x)(x). See, Venkataramanan et al. 2002, and E. J.Fordham, A. Sezginer, and L. D. Hall, Imaging multiexponentialrelaxation in the (y, log _(e) T ₁) plane, with application to clayltration in rock cores, J. Mag. Resonance A, 113:139-150, 1995. Next, itis also mathematically ill-conditioned. Often, regularization or priorinformation about the expected solution is incorporated into the problemformulation to make it better conditioned. However, the choice of theregularization functional as well as the weight given to the priorinformation is a well-known drawback of this transform due to itsnon-uniqueness. See, W. H. Press, S. A. Teukolsky, and W. T. Vetterling,Numerical Recipes in C, Cambridge University Press, 1992. Thus errors inthe estimation of T₁, T₂, or D distributions are propagated into errorsin the estimation of fluid properties.

SUMMARY

According to some embodiments, a method for estimating a property of amaterial is provided. A measurable change in the material is induced andmeasurements are made of an aspect of the material that tends toexponentially decay following the induction. The induction is repeatedat predetermined intervals, and the material property is estimated basedon measurements of the aspect. The intervals are designed so as toimprove the estimation of the property.

The induction can be a transmitted transverse magnetic pulse as with NMRmeasurements, according to some embodiments. According to someembodiments the induction is the induction of a pressure or temperaturechange in the material. According to some embodiments, the intervalspacing (e.g. in time, in frequency, etc.) is designed according to aformula determined by the property being estimated, and the spacing isnon-uniform in both linear and logarithmic scales. According to someembodiments the interval spacing in time follows a power law and allowsfor the material property to be directly computed from the measurements.

The method can be carried out as part of an oilfield service, and themeasurements can be made downhole. The material can be a fluid, a solidor a mixture. According to some embodiments, the material is one or morerock solids and the material properties being estimated arepetrophysical properties of the rock solids.

According to some embodiments, a system for estimating a property of amaterial is provided that includes an inducer that is adapted to inducea measurable change in the material, and one or more sensors adapted tomeasure an aspect of the material that tends to exponentially decayfollowing an induction. A controller is coupled to the inducer, andadapted and programmed to repeat the induction at predeterminedintervals. A processing system is adapted to estimate the materialproperty based on measurements. The induction intervals are designed soas to improve the estimation of the property.

Further features and advantages will become more readily apparent fromthe following detailed description when taken in conjunction with theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is further described in the detailed descriptionwhich follows, in reference to the noted plurality of drawings by way ofnon-limiting examples of exemplary embodiments, in which like referencenumerals represent similar parts throughout the several views of thedrawings, and wherein:

FIG. 1 is a plot showing the values of n and μ for different values of ωas determined by equation (8) herein;

FIG. 2 is a plot showing an example of a single exponential acquired atuniform time intervals in the time-domain;

FIG. 3 is a plot showing a non-uniform sampling, according to someembodiments;

FIG. 4 plots examples of non-uniform sampling in time domain that candirectly provide moments of relaxation or diffusion distribution,according to some embodiments;

FIG. 5 is a plot showing the ratio of the error-bar in the momentsobtained from uniform sampling to non-uniform sampling;

FIG. 6A shows an NMR tool being deployed in a wellbore, according toembodiments;

FIG. 6B is a schematic of a processing system used to control, recordand process data from an NMR tool; and

FIG. 7 illustrates another example of a wellsite system in which an NMRtool can be employed, according to some embodiments.

FIG. 8 represents magnetization data simulated with additive white noisewith σ_(ε)=0.01 at three different temperatures.

FIG. 9 represents T₂ distributions obtained for data shown in FIG. 1using inverse Laplace transform with regularization described inequation (6).

FIG. 10 represents moments obtained for data shown in FIG. 1 usingequations. (3)-(7)

FIG. 11 represents magnetization data simulated with additive whitenoise with σ_(s)=0.01 at three different temperatures.

FIG. 12 shows T₂ distributions obtained for data shown in FIG. 4 usinginverse Laplace transform with regularization described in equation (6).

FIG. 13 shows moments obtained for data shown in FIG. 4 using equations.(3)-(7).

FIG. 14 represents for given ω that the values of n and μ are unique anddetermined by eqn. (7).

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The following description provides exemplary embodiments only, and isnot intended to limit the scope, applicability, or configuration of thedisclosure. Rather, the following description of the exemplaryembodiments will provide those skilled in the art with an enablingdescription for implementing one or more exemplary embodiments. It beingunderstood that various changes may be made in the function andarrangement of elements without departing from the spirit and scope ofthe invention as set forth in the appended claims.

Specific details are given in the following description to provide athorough understanding of the embodiments. However, it will beunderstood by one of ordinary skill in the art that the embodiments maybe practiced without these specific details. For example, systems,processes, and other elements in the invention may be shown ascomponents in block diagram form in order not to obscure the embodimentsin unnecessary detail. In other instances, well-known processes,structures, and techniques may be shown without unnecessary detail inorder to avoid obscuring the embodiments. Further, like referencenumbers and designations in the various drawings indicated likeelements.

Also, it is noted that individual embodiments may be described as aprocess which is depicted as a flowchart, a flow diagram, a data flowdiagram, a structure diagram, or a block diagram. Although a flowchartmay describe the operations as a sequential process, many of theoperations can be performed in parallel or concurrently. In addition,the order of the operations may be re-arranged. A process may beterminated when its operations are completed, but could have additionalsteps not discussed or included in a figure. Furthermore, not alloperations in any particularly described process may occur in allembodiments. A process may correspond to a method, a function, aprocedure, a subroutine, a subprogram, etc. When a process correspondsto a function, its termination corresponds to a return of the functionto the calling function or the main function.

Furthermore, embodiments of the invention may be implemented, at leastin part, either manually or automatically. Manual or automaticimplementations may be executed, or at least assisted, through the useof machines, hardware, software, firmware, middleware, microcode,hardware description languages, or any combination thereof. Whenimplemented in software, firmware, middleware or microcode, the programcode or code segments to perform the necessary tasks may be stored in amachine readable medium. A processor(s) may perform the necessary tasks.

According to some embodiments, modified pulse sequences are providedthat allow for the direct measurements of moments of relaxation-time ordiffusion distributions. The data from these pulse sequences providesdirect estimates of fluid properties such as average chain length andviscosity of a hydrocarbon. In comparison to the uniformly-spaced pulsesequence, these pulse sequences are faster and have a lower error bar incomputing the fluid properties.

Co-pending U.S. Patent Application No. 61/242,218, filed Sep. 14, 2009(incorporated herein by reference and copied below) describes a methodto analyse NMR data that obviates the estimation of the T₁, T₂, or Ddistributions.

This invention relates to the field of estimation of sample properties.One preferred application of such scheme would be estimation of theproperties of a sample from exponentially decaying data.

Fitting exponentials to measured data is a well-known ill-conditionedproblem in science and engineering. It involves solving for non-negativeamplitude with time constant T₂ and amplitudes ƒ_(T) ₂ from the measuredmulti-exponential decay M(t) in the time-domain:

$\begin{matrix}{{M(t)} = {\int_{0}^{\infty}{e^{{- t}/T_{2}}{f_{T_{2}}\left( T_{2} \right)}{{dT}_{2}.}}}} & (1)\end{matrix}$

The time constants T₂ are often assumed to be a continuum. Without lossof generality, the corresponding non-negative amplitude ƒ_(T) ₂ isconsidered to be the probability density function of variable T₂.

Nuclear magnetic resonance (NMR) applications in biological systems,porous media or well-logging are further discussed hereinafter asexemplary applications comprising exponentially decaying data. However,the method of the invention relates to any type of exponentiallydecaying data.

Traditional well-logging using Nuclear Magnetic Resonance (NMR) isemployed to infer petrophysical properties of the formation as well asin-situ fluid properties. The multi-exponential time-decay of NMRmagnetization is characterized by relaxation times T₁ or T₂, andcorresponding amplitudes ƒ_(T) ₁ (T) or ƒ_(T) ₂ (T₂) (See R. L.Kleinberg, C. Straley, W. E. Kenyon, R. Akkurt, and S. A. Farooqui.Nuclear magnetic resonance of rocks: T1 vs. T2. SPE, 26470:553-563,1993). Often, the magnetization data has a wide range of relaxationtimes and ƒ_(T) ₁ (T₁) or ƒ_(T) ₂ (T₂) corresponds to a distribution ofthe underlying relaxation times.

These distributions provide information such as the range of pore-sizesin a rock and are empirically related to rock permeability. They arealso used to distinguish bound fluid versus free-fluid in a rock. NMRrelaxation measurements of bulk hydrocarbons are also used to infer oilviscosity. More recently, two-dimensional measurements, such asdiffusion-relaxation, have been used to provide information about fluidtyping and saturation (M. D. Hurlimann and L. Venkataramanan.Quantitative measurement of two dimensional distribution functions ofdiffusion and relaxation in grossly homogeneous fields. Journal MagneticResonance, 157:31-42, July 2002).

The core of the NMR interpretation relies on the inverse Laplacetransform of the magnetization data to reconstruct ƒ_(T) ₁ (T) or ƒ_(T)₂ (T₂). It is well-known that computation of this transform is amathematically ill-conditioned problem (L. Venkataramanan, Y.-Q. Song,and M. D. Hurlimann. Solving Fredholm integrals of the first kind withtensor product structure in 2 and 2.5 dimensions. IEEE Transactions onSignal Processing, 50:1017-1026, 2002). Often, regularization or priorinformation about the expected solution needs to be incorporated to makethe problem better conditioned. However, the choice of theregularization functional as well as the weight given to priorinformation is a well-known drawback of this transform due to itsnon-uniqueness.

In oilfield applications, the measured NMR magnetization data denoted byM(t) is a multiexponential decay, with time constant T₂ and amplitudesƒ_(T) ₂ .

The relaxation time T₂ is the characteristic time corresponding to lossof energy by protons in hydrocarbons or water present in pores of a rockor in the bulk fluid. There are a number of mechanisms for loss of thisenergy: frequent collision of molecules with the pore surface (surfacerelaxation), diffusion in a magnetic field (diffusion-inducedrelaxation) and interaction with other fluid molecules (bulkrelaxation).

As seen in eqn. (1), the measured data M(t) is a Laplace transform ofthe amplitudes ƒ_(T) ₂ (T₂). Traditionally, the inverse Laplacetransform is used to estimate the time constants T₂ and amplitudes ƒ_(T)₂ (T₂) from the measured data. The time constants T₂ are often assumedto be a continuum; thus the corresponding amplitudes ƒ_(T) ₂ (T₂) arereferred to as the T₂ distribution. The amplitude ƒ_(T) ₂ (T₂) at anygiven T₂ is proportional to the number of protons relaxing at that rate.Thus, the area under the T₂ distribution is interpreted as the porosityof the rock. Often, based on lithology, a threshold T₂ is chosen as thecut-off characteristic time separating fluid bound to the pore surfaceand fluid that is not bound to the pore surface and can flow moreeasily. For example, in sandstones, the area under the T₂ distributioncorresponding to relaxation times smaller than 33 msec has beenempirically related to bound fluid volume. The remaining area under theT₂ distribution corresponds to free fluid volume. The mean of thedistribution, <ln(T₂)> is empirically related to either rockpermeability and/or to hydrocarbon viscosity. The width of thedistribution, σ_(ln(T) ₂ ₎ provides a qualitative range of thedistribution of pore sizes in the rock. Similar answer products obtainedfrom either areas or moments are computed from two-dimensionaldiffusion-relaxation data or T₁-T₂ relaxation data (See Y.-Q. Song, L.Venkataramanan, M. D. Hurlimann, M. Flaum, P. Frulla, and C. Straley. T1T2 correlation spectra obtained using a fast two-dimensional Laplaceinversion. Journal Magnetic Resonance, 154:1-8, 2002, and D. Allen, S.Crary, and R. Freedman. How to use borehole Nuclear Magnetic Resonance.Oilfield Review, pages 34-57, 1997).

Eqn. (1) provides the forward model for T₂ relaxation. Other pulsesequences used to measure petrophysical or fluid properties such asdiffusion D or longitudinal relaxation T₁ result in mathematicalformulations similar to eqn. (1), with the variable T₂ being replaced byrelaxation T₁ or 1/D. Thus, analysis developed for solving eqn. (1) forƒ_(T) ₂ (T₂) can also be used to solve for distributions of longitudinalrelaxation ƒ_(T) ₁ (T₁) and diffusion f_(D)(D).

It is well known in the literature that the inverse Laplace transform isan ill-conditioned problem. Small changes in the measured data due tonoise can result in widely different T₂ distributions. In theory, thereare infinitely many solutions for ƒ_(T) ₂ (T₂). In the literature, thisproblem is typically solved by subjectively choosing the “smoothest”solution ƒ_(T) ₂ (T₂) that fits the data. This smooth T₂ distribution isoften estimated by minimization of a cost function Q with respect to theunderlying ƒ (E. J. Fordham, A. Sezginer, and L. D. Hall. Imagingmultiexponential relaxation in the (y, log_(e)T₁) plane, withapplication to clay filtration in rock cores. J. Mag. ResonanceA,113:139-150, 1995),Q=∥M−Kf∥ ²+α∥ƒ∥²,  (2)where M is the measured data, K is the matrix of the discretized kernele^(−t/T) ² and ƒ is the discretized version of the underlying densityfunction ƒ_(T) ₂ (T₂) in eqn. (1). The first term in the cost functionis the least squares error between the data and the fit from the modelin eqn. (1). The second term is referred to as regularization andincorporates smoothness in the relaxation amplitudes into the problemformulation. The parameter α denotes the compromise between the fit tothe data and an a priori expectation of the distribution.

The use of regularization and incorporation of smoothness into theproblem formulation is subjective. In eqn. (2), it is possible to definealternate cost functions based on entropy, the Backus-Gilbert method orother measures of difference between data and fit (W. H. Press, S. A.Teukolsky, and W. T. Vetterling. Numerical Recipes in C. CambridgeUniversity Press, 1992). In the next section, we describe an alternatemethod that obviates the need for the use of the inverse Laplacetransform and the computation of the T₂ distribution to estimateparameters such as

ln(T₂)

and σ_(ln(T) ₂ ₎

An alternate method is described that operates directly on the measuredmagnetization data and provides information about a sample parameterswithout the need for this ill conditioned inverse Laplace transform.

Consider a function of the form:G(ω)≡ln

T ₂ ^(ω)

  (3)where ω is a real number,

.

refers to the expectation operator and

T₂

refers to the expected value of the ω-th moment of T₂, defined as

$\begin{matrix}{{\left\langle T_{2}^{\omega} \right\rangle \equiv \frac{\int_{0}^{\infty}{T_{2}^{\omega}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}{\int_{0}^{\infty}{{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}} = {\frac{1}{\phi}{\int_{0}^{\infty}{T_{2}^{\omega}{F_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}}} & (4)\end{matrix}$where ϕ is the porosity and

ϕ = M(t = 0) = ∫₀^(∞)f_(T₂)(T₂)dT₂.By construction, G( ) is the second characteristic function of ln(T₂)(see Appendix A), with

$\begin{matrix}{{{{{\frac{dG}{d\;\omega}}_{\omega = 0} = \left\langle {\ln\left( T_{2} \right)} \right\rangle},\frac{d^{2}G}{d\;\omega^{2}}}}_{\omega = 0} = {\sigma_{l\;{n{(T_{2})}}}^{2}.}} & (5)\end{matrix}$

In eqn. (3), G(ω) can be computed from moments

T₂ ^(ω)

of relaxation time T₂ which, in turn, can be computed directly from themeasurement (see Appendix B),

$\begin{matrix}{{\left\langle T_{2}^{\omega} \right\rangle = {\frac{\left( {- 1} \right)^{n}}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d^{n}{M(t)}}{{dt}^{n}} \right\rbrack}{dt}}}}},{\omega = {\mu - n}}} & (6)\end{matrix}$where Γ( ) represents the Gamma function. The contribution of variable ωis in two parts: a real number μ and an integer n where the mathematicaloperator t^(μ-1) operates on the n-th derivative of the data. As shownin Appendix B, whenω>0:n=0,μ=ωω≤0:ω=μ−n,0<μ≤1,n=[−ω]+1.  (7)where [ω] refers to the integral part of the number ω. In eqn. (1), themoments

T₂ ^(ω)

are computed as a weighted linear combination of the measured data.This, in turn, allows easy computation of the function G( ) in eqn. (3)and thus, the parameters

ln(T₂)

and σ_(ln(T) ₂ ₎ in eqn. (5).

In another embodiment of the invention, the moments of relaxation timecan also be calculated in the Fourier domain (see Appendix C).

EXAMPLES OF APPLICATIONS

In addition to direct computation of parameters such as

ln(T₂)

and σ_(ln(T) ₂ ₎, the computation of moments has numerous applications.Some of these applications are listed below.

1 Hydrocarbon Properties

We demonstrate below that the computation of moments of relaxation timeis useful for example in the computation of hydrocarbon properties suchas average chain length, viscosity and parameters relating fluidcomposition to NMR data.

Scaling laws established on mixtures of alkanes provide a relationbetween fluid composition and relaxation time and/or diffusioncoefficients of the components of the mixture (See D. E. Freed. Scalinglaws for diffusion coefficients in mixtures of alkanes. Physical ReviewLetters, 94:067602, 20058, and D. E. Freed. Dependence on chain lengthof NMR relaxation times in mixtures of alkanes. Journal of ChemicalPhysics, 126:174502, 2007). A typical crude oil has a number ofcomponents, such as methane, ethane etc. Let N_(i) refer to the numberof carbon atoms or the chain length of the i-th component. Let N denotethe molar average chain length of the crude oil, defined as the harmonicmean of the chain lengths,

$\begin{matrix}{{\frac{1}{\overset{\_}{N}} \equiv {\sum\limits_{i = 1}^{\infty}\frac{f_{N}\left( N_{i} \right)}{N_{i}}}} = \left\langle N^{- 1} \right\rangle} & (8)\end{matrix}$where ƒ_(N) (N) denotes the mass fraction of component i in thehydrocarbon. According to the scaling law, the bulk relaxation timeT_(2,i) of component i in a crude oil follows a power law and dependsinversely on its chain length or number of carbon atoms, N_(i). Thus,the larger the molecule, the smaller the relaxation time T_(2,i). Therelaxation time also depends inversely on the average chain length N ofthe hydrocarbon, according toT _(2,i) =BN _(i) ^(−κ) N ^(−γ).  (9)

In eqn. (9), parameters B, γ, κ are constants at a given temperature andpressure. Similarly, if D refers to the diffusion coefficient ofcomponent iD ₁ =AN _(i) ^(−ν) N ^(−β)  (10)where A,ν,β are constants at a given temperature and pressure. Thescaling laws in eqns. (9) and (10) imply that measurements of diffusionor relaxation can be used to characterize fluid composition. The scalingtheory can also be used to derive an expression for fluid viscosity,η∝

D ⁻¹

η∝N ^(β)

N ^(ν)

  (11)

Traditional analysis of NMR relaxation data to compute either theaverage chain length or viscosity involves the following steps. First,one of the infinitely many solutions for ƒ_(T) ₂ (T₂) is computed fromeqn. (2). Next, the chain length distribution ƒ_(N) (N_(i)) is computedusing eqn. (9). Next, fluid properties such as the average chain lengthand viscosity are computed using eqn. (11). Errors in the estimation ofƒ_(T) ₂ (T₂) from the ill-conditioned nature of the inverse Laplacetransform are propagated into the estimation of fluid properties.

We show below that the computation of moments of either the relaxationtimes or diffusion constants allows direct estimation of the fluidproperties. Taking the expectation operator on either side of Eqns. (9)and (10) and re-writing it, the average chain length from either therelaxation-time or diffusion data is

$\begin{matrix}{\overset{\_}{N} = \left\lbrack {B^{{- 1}/\kappa}\left\langle T_{2}^{1/\kappa} \right\rangle} \right\rbrack^{- \frac{\kappa}{\gamma + \kappa}}} & (12) \\{\overset{\_}{N} = {\left\lbrack {A^{{- 1}/v}\left\langle D^{1/v} \right\rangle} \right\rbrack^{\frac{v}{v + \beta}}.}} & (13)\end{matrix}$

From eqn. (11), the viscosity of the hydrocarbon can be computeddirectly from the moments of relaxation time,

$\begin{matrix}{\eta \propto {\left\lbrack \left\langle T_{2}^{\frac{1}{\kappa}} \right\rangle \right\rbrack^{\frac{{\beta\mspace{11mu}\kappa} - {\gamma\; v}}{\gamma + \kappa}}{\left\langle T_{2}^{- \frac{v}{\kappa}} \right\rangle.}}} & (14)\end{matrix}$

Thus, using eqns. (12) and (14) fluid properties such as average chainlength and the fluid viscosity of a hydrocarbon can be computed directlyfrom the moments of the relaxation time, which, in turn, can be directlycomputed from the data using eqns. (6) and (7). Similarly, eqns. (11)and (13) enable direct computation of fluid properties from diffusiondata. In addition to characterizing the hydrocarbon, the average chainlength can be a useful parameter to help classify hydrocarbons (M. D.Hurlimann, D. E. Freed, L. J. Zielinski, Y. Q. Song, G. Leu, C. Straley,C. Cao Minh, and A. Boyd. Hydrocarbon composition from NMR diffusion andrelaxation data. SPWLA 49th logging symposium, 2008). The temperatureindependence of this parameter can also be used to identify changes inhydrocarbon properties with depth.

The parameters B,γ,κ in eqn. (9) and A,ν,β in eqn. (10) relaterelaxation times or diffusion coefficients to fluid composition in amixture. Given NMR and gas-chromatography (GC) data, we show below thatthe computation of G(ω) provides a direct method to estimate theseparameters in a crude oil. Taking the variance on either side of eqn.(9) yieldsσ_(ln(T) ₂ ₎ ²=κ²σ_(ln(N)) ²,  (15)where σ_(ln(T) ₂ ₎ ² can be computed from NMR data using eqn. (5) and aσ_(ln(N)) ² can be obtained from GC data which provides fluidcomposition. Thus, given NMR and GC data, κ for any crude oil isobtained as the ratio of the width or standard deviation of the T₂distribution to the width or standard deviation of the chain lengthdistribution in the log-domain. The parameters B and γ can be computedin a similar manner. Taking the expectation on either side of eqn. (9)yields

ln(T ₂)

=ln(B)−γ ln( N )−κ

ln(N)

.  (16)

The term

ln(T₂)

can be computed from NMR data using eqn. (5). Parameter κ is found fromeqn. (15). Parameters N and

ln(N)

are computed from GC data. Thus, eqn. (16) is linear in parameters ln Band γ. Given data from two or more hydrocarbons with similar properties,these parameters can be estimated. When the data are measured as afunction of temperature and pressure, parameters κ, B and γ can becomputed as a function of temperature and pressure. For example, frommeasured NMR and GC data, the value of κ on various crude oils in thelab obtained using eqn. (15) are shown in Table 1. The error-bar in κ isaround 0.05.

TABLE I Estimated κ values from different crude oils Sample 30 C. 60 C.90 C. 120 C. 1 1.46 1.38 1.29 1.22 2 1.51 1.46 1.41 N/A 3 1.37 1.29 1.221.16 4 1.54 1.39 1.34 1.21

For a given sample, the magnitude of κ seemingly decreases withincreasing temperature. The value of κ may also vary slightly from onesample to another. This may have an important implication in thecomputation of fluid composition from measured NMR data. Similaranalysis holds good for estimated parameters A, ν and β in eqn. (10).Given both gas-chromatography and NMR diffusion data,σ_(ln(D)) ²=ν²σ_(ln(N)) ²  (17)

ln(D)

=ln(A)−β ln( N )−ν

ln(N)

.  (18)

Thus, using eqns. (17) and (18), NMR diffusion and GC data can be usedto the estimate the parameters A, ν and β and study their variation withtemperature and pressure.

2 Identification of Fluid Phase Change

An abrupt change in fluid parameters with temperature and pressure canbe an indicator of fluid phase change such as wax or asphalteneprecipitation. For example, Table 2 shows the value of estimated κ withtemperature for a particular lab sample. It is seen that above 30 C thevalue of κ is around 1.4±0.1. Below 20 C, κ increases abruptly to1.7±0.1, indicating that the molecular dynamics has slowed, signaling aphase transition. In this lab example, wax was found to haveprecipitated in the hydrocarbon between 20 and 30 C. Thus monitoring thevalue of κ (or other moments of relaxation time or diffusion such as themean chain length) as a function of temperature or pressure can helpidentify abrupt changes in fluid properties indicating phasetransitions.

TABLE 2 Estimated κ values from one crude oil at different temperaturesindicates a phase transition between 20 and 30 C. Temperature (C.) κ 51.67 ± 0.1 10 1.67 ± 0.1 15 1.61 ± 0.1 20 1.75 ± 0.1 30  1.42 ± 0.01 501.30 ± 0.13 Computation of Error-Bars on Answer Products

The traditional analysis for computing T₂ distributions is intrinsicallynon-linear due to the non-negativity constraint on the distribution [3].Thus, computation of error-bars has been a longstanding issue. Since thenew analysis described in the previous section provides a linear andanalytical expression for the parameters of interest, the error-bars areeasily computed analytically. For example, when ω>0, using eqns. (6) and(7), we have

$\begin{matrix}{\sigma_{\langle T_{2}^{\omega}\rangle} = {\frac{\sigma_{ò}}{{\Gamma(\omega)}\phi}t_{E}^{\omega}\sqrt{\sum\limits_{n = 1}^{N_{e}}n^{{{2\omega} - 2}\;}}}} & (19)\end{matrix}$where σ_(ε) is the standard deviation of additive white noise in themeasurement, N_(e) is the number of echoes in the measured T₂ relaxationdata and t_(E) is the echo spacing. Eqn. (19) provides an insight intothe desired signal to noise ratio in the data required to produce adesired uncertainty in the computation of moments. In turn, this can beconverted into uncertainty in the desired fluid and petrophysicalparameters.

Computation of fluid properties along with associated uncertainty fromNMR data can play a significant role in determining reservoirconnectivity. Abrupt changes in fluid properties, such as fluidcomposition or viscosity between zones in a single well or betweenneighboring wells can, under some circumstances, be a marker forreservoir compartmentalization (See H. Elshahawi, L. Venkataramanan, D.McKinney, M. Flannery, O. C. Mullins, and M. Hashem. Combiningcontinuous fluid typing, wireline formation testers, and geochemicalmeasurements for an improved understanding of reservoir architecture.SPE Reservoir Evaluation and Engineering, pages 27-40, 2008). However,the abrupt differences in fluid properties can be identified reliablyonly when the significance of uncertainties is taken into account.Therefore, computation of fluid properties from NMR along withassociated uncertainties can help determine if fluids in neighboringzones are statistically similar or different. If the two fluids inadjacent layers are statistically different, there may be a permeabilitybarrier between the zones that divides the layer into separateflow-units.

Simulations

In the first simulation, NMR data were simulated from gas-chromatographyusing eqn. (9) with a constant value of κ=1.24 at three temperatures:10, 30 and 50 C. FIG. 8 shows one realization of data with additivewhite Gaussian noise with σ_(ε)=0.01.

The simulation model parameters are shown in columns 1-3 in Table 3. Thesimulated magnetization data are analyzed using two methods: thetraditional inverse Laplace transform described in [6] and the newanalysis. The resultant T₂ distributions obtained from the traditionalmethod are shown in FIG. 9. Column 4 in Table 3 provides the results ofthe analysis using this method. The error-bars were obtained usingMonte-Carlo analysis with a number of simulated data sets with differentnoise realizations. Comparing column 4 with column 3, we see that theparameters <ln(T₂)> and κ are biased, and associated with a largevariance. This bias results from estimated T₂ distributions that havebeen artificially smoothened by regularization. Artificial peaks atsmaller and larger relaxation times also contribute significantly to thebias and variance. FIG. 10 shows the moments obtained by the newanalysis. Column 5 in Table 3 provides the quantitative results on thedata. By comparing columns 3 and 5, it is seen that the new analysisperforms reasonably well. The parameters are estimated within atolerable bias and variance.

TABLE 3 Simulation 1: Comparison of results indicates that the newanalysis performs better than the traditional approach. Model True Inv.Laplace Parameters Temp(C.) Values Transform New Analysis

 log₁₀ (T₂) 

10 −0.036 −0.04 ± 0.01  −0.037 ± 0.002  30 0.032 0.023 ± 0.01  0.032 ±0.002 50 0.062 0.053 ± 0.01  0.061 ± 0.002 σ_(log10(T) ₂ ₎ ² 10 0.2 0.25± 0.06 0.198 ± 0.003 30 0.2 0.25 ± 0.06 0.199 ± 0.003 50 0.2 0.25 ± 0.060.200 ± 0.003 κ 10 1.24 1.54 ± 0.35 1.23 ± 0.02 30 1.24 1.55 ± 0.37 1.24± 0.02 50 1.24 1.56 ± 0.38 1.24 ± 0.03

In a second simulation, data were simulated from the same gaschromatography as the previous example of FIG. 8 with the exception thatκ decreases with increasing temperature as shown in column 3 in Table 4.The data are corrupted with white Gaussian noise and shown in FIG. 11.The results of the analysis using the traditional inverse Laplacetransform are shown in FIG. 12 and column 4 in Table 4. The results ofthe new analysis are shown in FIG. 13 and column 5. As before, theerror-bars were obtained from Monte-Carlo analysis with a number ofsimulated data sets with different noise realizations. Comparing columns3 and 5, it is seen that the new analysis performs better than theinverse Laplace transform; all the parameters have a smaller bias andvariance. Further, the trend in κ with temperature is also estimatedreasonably well.

TABLE 4 Simulation 2: Comparison of results indicates that the newanalysis performs better than the traditional approach. Model True Inv.Laplace Parameters Temp(C.) Values Transform New Analysis

 log₁₀ (T₂) 

10 −0.040 −0.05 ± 0.01  −0.040 ± 0.002  30 0.046 0.04 ± 0.01 0.040 ±0.002 50 0.057 0.052 ± 0.01  0.056 ± 0.002 σ_(log10(T) ₂ ₎ ² 10 0.210.26 ± 0.06  0.21 ± 0.003 30 0.18 0.22 ± 0.06  0.18 ± 0.004 50 0.11 0.16± 0.06  0.12 ± 0.005 κ 10 1.3 1.58 ± 0.34 1.28 ± 0.02 30 1.1 1.43 ± 0.391.11 ± 0.02 50 0.7 1.50 ± 0.38 0.75 ± 0.03

According to embodiment of the invention, moments of the relaxation timeand/or diffusion can be directly computed from the measuredmagnetization data. In turn, these moments are directly related toproperties of a sample like the petrophysical parameters such as

ln(T₂)

and σ_(ln(T) ₂ ₎. The moments can also directly provide fluid parameterssuch as the mean chain length N and viscosity η. Given NMR and GasChromatography data, other parameters that relate fluid composition torelaxation or diffusion data can also be estimated. The error-bars inthe estimated parameters can also be readily computed. In addition,monitoring of some of these fluid parameters with changes in temperatureand pressure can be an indicator of phase transitions in thehydrocarbon.

In this section, we show that the function G(ω) defined in eqn. (3) isthe second characteristic function of T₂. From eqn. (4), we can see thatthe function G(ω)=0 when ω=0.

$\begin{matrix}{{G\left( {\omega = 0} \right)} = {{{\ln\left\langle T_{2}^{\omega = 0} \right\rangle} \equiv {\ln\left\lbrack \frac{\int_{0}^{\infty}{T_{2}^{\omega = 0}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}{\int_{0}^{\infty}{{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}} \right\rbrack}} = 0}} & (20)\end{matrix}$

Next, eqn. (4) can be re-written as,

$\begin{matrix}{{\left\langle T_{2}^{\omega} \right\rangle \equiv \frac{\int_{0}^{\infty}{e^{{\omega ln}{(T_{2})}}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}{\int_{0}^{\infty}{{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}} = \frac{\int_{- \infty}^{\infty}{e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}{\int_{- \infty}^{\infty}{{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}} & (21)\end{matrix}$where ƒ_(ln(T) ₂ ₎ refers to the underlying density function of ln(T₂).Thus

$\begin{matrix}{\frac{{dG}(\omega)}{d\;\omega} = {\frac{d\;\ln\left\langle T_{2}^{\omega} \right\rangle}{d\;\omega} = \frac{\int_{- \infty}^{\infty}{{\ln\left( T_{2} \right)}e^{\omega\;\ln\;{(T_{2})}}{f_{\ln\;{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}{\int_{- \infty}^{\infty}{e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}}} & (22)\end{matrix}$and

$\begin{matrix}{\frac{d^{2}{G(\omega)}}{d\;\omega^{2}} = {\frac{d^{2}\ln\left\langle T_{2}^{\omega} \right\rangle}{d\;\omega^{2}} = {\frac{\int_{- \infty}^{\infty}{\left\lbrack {\ln\left( T_{2} \right)} \right\rbrack^{2}e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}{\int_{- \infty}^{\infty}{e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}} - \left( \frac{\int_{- \infty}^{\infty}{{\ln\left( T_{2} \right)}e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}}{\int_{- \infty}^{\infty}{e^{\omega\;{\ln{(T_{2})}}}{f_{\ln{(T_{2})}}\left( {\ln\left( T_{2} \right)} \right)}d\;{\ln\left( T_{2} \right)}}} \right)^{2}}}} & (23)\end{matrix}$

Therefore,

$\begin{matrix}{{\left. \frac{d\; G}{d\;\omega} \right|_{\omega = 0} = \left\langle {\ln\left( T_{2} \right)} \right\rangle},{\left. \frac{d^{2}G}{d\;\omega^{2}} \right|_{\omega = 0} = \sigma_{\ln{(T_{2})}}^{2}}} & (24)\end{matrix}$Similarly, higher-order moments or cumulants of ln(T₂) are obtained bytaking the higher-order derivatives of G( ) with respect to ω.

Appendix B

In this section, we show that moments of ln(T₂) defined in eqn. (4) canbe computed directly from the measurement, from the definition of Γ(.)function. Consider the following two cases.

Case 1: When ω>0. The Γ(.) function is a generalization of the factorialfunction and is defined as,

$\begin{matrix}{{\Gamma(\omega)} \equiv {\int_{0}^{\infty}{t^{\omega - 1}e^{- t}{dt}}}} & (25)\end{matrix}$

Therefore, when the operator t^(ω-1) operates on a single exponentiale^(−at) where α>0,

$\begin{matrix}{{\int_{0}^{\infty}{t^{\omega - 1}e^{{- a}\; t}{dt}}} = \frac{\Gamma(\omega)}{a^{\omega}}} & (26)\end{matrix}$

When the operator t^(ω-1) operates on a sum of exponentials given byeqn. (1), then,

$\begin{matrix}\begin{matrix}{{\int_{0}^{\infty}{t^{\omega - 1}{M(t)}{dt}}} = {\int_{0}^{\infty}{t^{\omega - 1}{\int_{0}^{\infty}{e^{{- t}/T_{2}}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}{dt}}}}}} \\{= {\int_{0}^{\infty}{\left\lbrack {\int_{0}^{\infty}{t^{\omega - 1}e^{{- t}/T_{2}}{dt}}} \right\rbrack{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}} \\{= {{\Gamma(\omega)}{\int_{0}^{\infty}{T_{2}^{\omega}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}}}\end{matrix} & (27)\end{matrix}$

From eqns. (4) and (27),

$\begin{matrix}{\left\langle T_{2}^{\omega} \right\rangle = {\frac{1}{{\Gamma(\omega)}\phi}{\int_{0}^{\infty}{t^{\omega - 1}{M(t)}{dt}}}}} & (28)\end{matrix}$

Case 2: When ω≤0.

Consider the integrand

$t^{\mu - 1}\frac{d^{n}M}{{dt}^{n}}$where n≥0 is an integer and 0<μ≤1. Let ω=μ−n. From eqn. (1),

$\begin{matrix}\begin{matrix}{{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d^{n}{M(t)}}{{dt}^{n}} \right\rbrack}{dt}}} = {\left( {- 1} \right)^{n}{\int_{0}^{\infty}{\left\lbrack {\int_{0}^{\infty}{\frac{t^{\mu - 1}}{T_{2}^{n}}e^{{- t}/T_{2}}}} \right\rbrack{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}}} \\{= {\left( {- 1} \right)^{n}{\Gamma(\mu)}\left\langle T_{2}^{\mu - n} \right\rangle\phi}} \\{= {\left( {- 1} \right)^{n}{\Gamma(\mu)}\left\langle T_{2}^{\omega} \right\rangle\phi}}\end{matrix} & (29)\end{matrix}$

Thus

$\begin{matrix}{\left\langle T_{2}^{\omega} \right\rangle = {\frac{\left( {- 1} \right)^{n}}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d^{n}{M(t)}}{{dt}^{n}} \right\rbrack}{dt}}}}} & (30)\end{matrix}$

Given a value of ω≤0, there are a number of different ways to compute μand n such that the equation ω=μ−n where n is an integer is valid. Oneof these methods is to choose n and μ such thatn=[−ω]+1 and 0<μ≤1  (31)where [ω] refers to the integral part of real number ω. Usingintegration by parts, it can be shown that all the other methods ofchoosing n and μ collapse to eqn. (31). FIG. 14 shows the variation of nand μ for the different values of ω.

In another aspect of the invention, porosity can also be estimateddirectly from the measured data. Let R(ω)=ϕ

T₂ ^(ω)

. From eqn. (6),

$\begin{matrix}{{{R(\omega)} = {\frac{\left( {- 1} \right)^{n}}{\Gamma(\mu)}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d^{n}{M(t)}}{{dt}^{n}} \right\rbrack}{dt}}}}},{{\mu - n} = \omega}} & (32)\end{matrix}$

where the relation between n, μ and ω is given in eqn. (7). Bydefinition in eqn. (4),

T₂ ^(ω=0)

=1. Therefore, the porosity can be computed as ϕ=R(ω=0), whichcorresponds to μ=n=1.

Appendix C

In this section, we show another embodiment of the invention to computethe moments of relaxation time and/or diffusion. This embodiment takesadvantage of the fact that the Fourier transform of a Gaussian is aGaussian with a width in the Fourier domain that is inverselyproportional to the width in the time domain. The moments can becomputed directly from the data,

$\begin{matrix}{\left\langle T_{2}^{\omega} \right\rangle = \left\{ \begin{matrix}{\frac{1}{{\Gamma(\omega)}\phi}{\int_{0}^{\infty}{t^{\omega - 1}{M(t)}{dt}}}} & {{{if}\mspace{14mu}\omega} > 0} \\{\frac{2\pi^{({0.5 - {2\;\omega}})}}{{\Gamma\left( {0.5 - \omega} \right)}\phi}{\int_{0}^{\infty}{k^{{- 2}\;\omega}{{\overset{\sim}{M}}_{y}(k)}{dk}}}} & {{{if}\mspace{14mu}\omega} < {0.5.}}\end{matrix} \right.} & (33)\end{matrix}$where {tilde over (M)}_(y)(k) is the Fourier transform of the measureddata in √{square root over (t)} domain,

$\begin{matrix}{{{{\overset{\sim}{M}}_{y}(k)} = {\int_{0}^{\infty}{{M_{y}(y)}e^{{- 2}\;\pi\;{iky}}{dy}}}},{{{where}\mspace{14mu}{M_{y}(y)}} = {M\left( {t = y^{2}} \right)}}} & (34)\end{matrix}$The proof of eqns. (33) and (34) is given below.Case 1: When ω>0: The proof of the first part of eqn. (33) is given ineqns. (25)-(28) in

Appendix B.

Case 2: When ω<0.5: The proof of this section depends on the followingeqns. (35)-(38). Using eqn. (25), it is seen that∫₀ ^(∞) t ^(μ-1) e ^(−t) ² dt=½Γ(μ/2)wenμ>0  (35)Therefore, for any α>0 and μ>0,

$\begin{matrix}{{\int_{0}^{\infty}{t^{\mu - 1}e^{- {at}^{2}}{dt}}} = {\frac{1}{2a^{\mu/2}}{\Gamma\left( {\mu\text{/}2} \right)}}} & (36)\end{matrix}$

Let {tilde over (F)}(k) denote the Fourier transform of a function ƒ(t),

$\begin{matrix}{{\overset{\sim}{F}(k)} = {\int_{- \infty}^{\infty}{{f(t)}e^{{- 2}\pi\; i\;{kt}}{dt}}}} & (37)\end{matrix}$

The Fourier transform of a Gaussian is a Guassian with a width in thefrequency domain that is inversely proportional to the width in the timedomain. Thus when

$\begin{matrix}{{{f(t)} = e^{- {ct}^{2}}},{{{then}\mspace{14mu}{\overset{\sim}{F}(k)}} = {\sqrt{\frac{\pi}{c}}e^{\frac{{- \pi^{2}}k^{2}}{c}}}}} & (38)\end{matrix}$Let M_(y)(y)=M(t=y²) where M(t) is the measured data given in eqn. (1),

$\begin{matrix}{{M_{y}(y)} = {\int_{0}^{\infty}{e^{- \frac{y^{2}}{T_{2}}}{f_{T_{2}}\left( T_{2} \right)}{{dT}_{2}.}}}} & (39)\end{matrix}$

Let {tilde over (M)}_(y)(k) denote the Fourier transform of M_(y),

$\begin{matrix}{{{\overset{\sim}{M}}_{y}(k)} = {\int_{- \infty}^{\infty}{{M_{y}(y)}e^{{- 2}\pi\;{iky}}{dy}}}} & (40)\end{matrix}$

From eqns. (39) and (40),

$\begin{matrix}{{{\overset{\sim}{M}}_{y}(k)} = {\int_{0}^{\infty}{\left\lbrack {\int_{- \infty}^{\infty}{e^{- \frac{y^{2}}{T_{2}}}e^{{- 2}\;\pi\;{iky}}{dy}}} \right\rbrack{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}} & (41)\end{matrix}$

Using eqns. (38) and (41)

$\begin{matrix}{{{\overset{\sim}{M}}_{y}(k)} = {\sqrt{\pi}{\int_{0}^{\infty}{\sqrt{T_{2}}e^{{- \pi^{2}}k^{2}T_{2}}{f_{T_{2}}\left( T_{2} \right)}{dT}_{2}}}}} & (42)\end{matrix}$

Now,

$\begin{matrix}{{\int_{0}^{\infty}{k^{\mu - 1}{{\overset{\sim}{M}}_{y}(k)}{dk}}} = {\sqrt{\pi}{\int_{0}^{\infty}{\left\lbrack {\int_{0}^{\infty}{k^{\mu - 1}e^{{- \pi^{2}}k^{2}T_{2}}{dk}}} \right\rbrack\sqrt{T_{2}}{f\left( T_{2} \right)}{dT}_{2}}}}} & (43)\end{matrix}$

Using eqns. (36) with α=π²T₂ and (43)

$\begin{matrix}{{\int_{0}^{\infty}{k^{\mu - 1}{M_{y}(k)}{dk}}} = {{\frac{(\pi)^{0.5 - \mu}{\Gamma\left( {\mu\text{/}2} \right)}}{2}\left\langle T_{2}^{0.5{({1 - \mu})}} \right\rangle\phi\mspace{14mu}{when}\mspace{14mu}\mu} > 0}} & (44)\end{matrix}$

Let 0.5(1−μ)=ω. When μ>0, ω<0.5. Therefore, for ω<0.5

$\begin{matrix}{\left\langle T_{2}^{\omega} \right\rangle = {\frac{2\pi^{({0.5 - {2\omega}})}}{{\Gamma\left( {0.5 - \omega} \right)}\phi}{\int_{0}^{\infty}{k^{{- 2}\omega}{{\overset{\sim}{M}}_{y}(k)}{dk}}}}} & (45)\end{matrix}$where {tilde over (M)}_(y) is the Fourier transform of the magnetizationdata M_(y) in eqn. (40) and M_(y) is the magnetization data in the√{square root over (t)} domain.

In another aspect of the invention, porosity can also be estimateddirectly from the measured data in the Fourier domain. Let R(ω)=ϕ

T₂ ^(ω)

. From eqn. (45),

$\begin{matrix}{{R\left( {\omega = 0} \right)} = {\phi = {\frac{2\sqrt{\pi}}{\Gamma(0.5)}{\int_{0}^{\infty}{{{\overset{\sim}{M}}_{y}(k)}{dk}}}}}} & (46)\end{matrix}$

A new mathematical formulation is described where the moments of arandom variable can be directly computed from the measurement. Themethod has two distinct advantages: the computation of moments is linearand does not require prior knowledge or information of the expecteddistribution. A moment of a random variable x with density functionƒ_(x)(x) is defined as,

$\begin{matrix}{\left\langle x^{\omega} \right\rangle \equiv \frac{\int_{0}^{\infty}{x^{\omega}{f_{x}(x)}{dx}}}{\int_{0}^{\infty}{{f_{x}(x)}{dx}}}} & (6)\end{matrix}$

For NMR applications, the variable x corresponds to relaxation times T₁,T₂, or diffusion D. The moments of relaxation times or diffusion can bedirectly computed from the measured magnetization data,

${\left\langle x^{\omega} \right\rangle = {\frac{\left( {- 1} \right)^{n}}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d^{n}{M(t)}}{{dt}^{n}} \right\rbrack}{dt}}}}},{\omega = {\mu - n}}$where Γ( ) represents the Gamma function and

ϕ = M(t = 0) = ∫₀^(∞)f_(x)(x)dx.The contribution of variable ω is in two parts: a real number μ and aninteger n where the mathematical operator t^(μ-1) operates on the n-thderivative of the data. Given a value of ω, the values of n and μ areunique and given as,ω>0: η=0,μ=ωω≤0:ω=μ−n,0<μ≤1,n=[−ω]+1  (8)where [ω] refers to the integral part of real number ω. FIG. 1 is a plot110 showing the values of n and μ for different values of ω asdetermined by above these equations.

According to some embodiments, the computation of moments of relaxationtimes T₁, T₂, or diffusion D, in conjunction with the scaling laws, canlead to greatly improved alternative pulse sequences to specificallytarget or estimate only certain moments of relaxation time/diffusion. Inco-pending U.S. Patent Application Publication No. 2008-0314582(incorporated herein by reference), “targeted” measurements aredescribed as measurements that target a certain specified parameter at aparticular location in the reservoir space. These targeted measurementsare achieved by both hardware and software modifications.

According to some embodiments, the concept of targeted measurements isapplied to NMR data and alternate pulse sequences are described todirectly estimate the moments of relaxation time or diffusioncoefficients. According to some embodiments, alternate pulse sequencesare designed by taking advantage of the specific exponential-kernelstructure in eqn. (1) as well as the computation of moments using eqn.(7). The magnetization data from the new pulse sequences can be directlyinterpreted to obtain fluid properties.

According to some embodiments, further detail will now be provided formodified pulse sequences for measuring NMR data. The non-uniformlyspaced pulse sequences are targeted to specifically measure one of themoments of the relaxation time T₁, T₂, or D distribution. Forsimplicity, consider that the random variable in eqn. (1) is T₁ and aspecific ω-th moment of T₁ is desired, where ω>0. This corresponds tothe right-hand side of the ω axis in FIG. 1. In this case, from eqns.(7) and (8),

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{1}{{\Gamma(\omega)}\phi}{\int_{0}^{\infty}{t^{\omega - 1}{M(t)}{dt}}}}} & (9)\end{matrix}$We introduce a new variable whereτ=t ^(ω)  (10)Therefore, eqn. (9) can be re-written as

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{1}{{\Gamma\left( {\omega + 1} \right)}\phi}{\int_{0}^{\infty}{{M_{\tau}(\tau)}d\;\tau}}}} & (11)\end{matrix}$where M_(τ)(τ) is the magnetization data in the τ domain. Eqn. (11)indicates that if the data are acquired such that they are uniformlyspaced in the τ domain, then the sum of the measured data directlyprovides the ω-th moment. From eqn. (10), uniform spacing in the τdomain implies non-uniform spacing in the time (t) domain.

We illustrate this with an example. FIG. 2 is a plot showing an exampleof a single exponential acquired at uniform time intervals in thetime-domain. In particular, curve 210 plots a single exponential decayof the form M(t)=e^(−t/T) ¹ in the time-domain as shown in where thedata are uniformly sampled in time. In this example, suppose that the0.5-th moment of T₁ is required. FIG. 3 is a plot showing a non-uniformsampling, according to some embodiments. In particular, using eqn. (10),the data shown by curve 310 are sampled uniformly in the √{square rootover (t)} domain. Curve 310 is an example of a single exponential in thetime-domain into a Guassian in the τ domain with M(t)=e^(−τ) ² ^(/T) ¹ .Using eqn. (11), the sum of measured echoes of the Gaussian directlyprovides

T₁ ^(0.5)

. Thus, if the ω-th moment is desired, then the data need to be samplednon-uniformly in time according to eqn. (10). FIG. 4 plots examples ofnon-uniform sampling in time domain that can directly provide moments ofrelaxation or diffusion distribution, according to some embodiments. Inparticular, plots 410, 412, 414 and 416 represent sampling for co valuesof 1, 0.75, 0.5, and 0.25 respectively. When 0<ω<1, this leads to finersampling at the initial part of the data and coarser sampling at thetail end of the data.

The concept of non-uniform sampling is also valid when ω<0. Consider thecase when −1<ω≤0. In this case, eqn. (7) reduces to

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {{\frac{- 1}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{{dM}(t)}{dt} \right\rbrack}{dt}\mspace{14mu}{where}\mspace{14mu}\mu}}} = {{1 + \omega} > 0}}} & (12)\end{matrix}$Eqn. (12) can be re-written as

$\begin{matrix}{{\left\langle T_{1}^{\omega} \right\rangle = {\frac{- 1}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 1}\left\lbrack \frac{d\left\lbrack {{M(t)} - {M(0)}} \right\rbrack}{dt} \right\rbrack}{dt}}}}}\;} & (13)\end{matrix}$Using integration by parts,

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {{\frac{- 1}{{\Gamma(\mu)}\phi}\left\lbrack {t^{\mu - 1}\left\{ {{M(t)} - {M(0)}} \right\}} \right\rbrack}_{t = 0}^{t = \infty} - {\frac{\left( {\mu - 1} \right)}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 2}\left\lbrack {{M(t)} - {M(0)}} \right\rbrack}{dt}}}}}} & (14) \\{\mspace{40mu}{= {{\frac{M(0)}{\phi}\delta_{\mu\; 1}} + {\frac{\left( {\mu - 1} \right)}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{{t^{\mu - 2}\left\lbrack {{M(t)} - {M(0)}} \right\rbrack}{dt}}}}}}} & (15)\end{matrix}$where δ_(μ1)=1 only if μ=1 and is zero otherwise. When μ=1(corresponding to ω=0), the first term in eqn. (15) is 1 (since M(0)=ϕ)and the second term is zero. This is consistent with the definition ofthe zero-th moment in eqn. (6) with

T₁ ^(ω=0)

=1.

When 0<μ<1, the first term in eqn. (15) is zero. The second term isintegrable at t=0 since [M(t)−M(0)]˜t as t→0. In this case, we introducea new variable τ whereτ=t ^(μ=1)  (16)Here, τ→∞ when t→0, and τ→0 when t→∞. Therefore, eqn. (15) becomes

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{- 1}{{\Gamma(\mu)}\phi}{\int_{0}^{\infty}{\left\lbrack {{M_{\tau}({`\tau})} - {M(0)}} \right\rbrack d\;\tau}}}} & (17)\end{matrix}$Similar to eqn. (11), the area under the integrand [M_(T)−M(0)] can beused to directly estimate the ω-th moment. As shown herein below, theconcept of the pulse sequence for ω>0 and −1<ω≤0 is valid for otherranges of ω as well. For example, when −2<ω≤−1, pulse sequences thatdirectly measure the ω-th moment are obtained by Taylor series expansionof the data with second-order terms in eqn. (13). Higher-order terms inthe Taylor series lead to further negative values of ω. As shown hereinbelow, when ω is not a negative integer,

T₁ ¹⁰⁷

can be obtained by sampling uniformly in the τ domain, whereτ=t ^(ω)  (18)When ω is a negative integer, from eqns. (7) and (8), the ω-th moment isobtained by the |ω|-th derivative of the data,

$\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{- 1^{\omega }}{\phi}\frac{d^{\omega }{M\left( {t = 0} \right)}}{{dt}^{\omega }}}} & (19)\end{matrix}$

In addition to computing the moments directly, this modified pulsesequence also enables fluid properties such as the average chain lengthand viscosity to be directly estimated in a simple manner. Re-writingeqns. (3) and (4), the average chain length and viscosity can becomputed from moments of T₁ or T₂ relaxation time or diffusion, as

$\begin{matrix}{\overset{\_}{N} = \left\lbrack {B^{{- 1}/\kappa}\left\langle T_{1}^{1/\kappa} \right\rangle} \right\rbrack^{- \frac{\kappa}{\gamma + \kappa}}} & (20) \\{\overset{\_}{N} = \left\lbrack {A^{{- 1}/v}\left\langle T_{1}^{1/v} \right\rangle} \right\rbrack^{- \frac{v}{v + \beta}}} & (21) \\\left. {\eta \propto \left\langle D^{- 1} \right\rangle}\Rightarrow{\eta \propto {\left\lbrack \left\langle T_{1}^{\frac{1}{\kappa}} \right\rangle \right\rbrack^{\frac{{\beta\;\kappa} - {\gamma\; v}}{\gamma + \kappa}}\left\langle T_{1}^{- \frac{v}{\kappa}} \right\rangle}} \right. & (22)\end{matrix}$According to eqns. (20) and (21), the average chain length can beestimated from a modified pulse sequence to directly estimate the(1/κ)-th moment of relaxation time or (1/ν)-th moment of diffusion.Similarly, viscosity can be estimated from two modified pulse sequencesto directly estimate the (1/κ)-th and (1/ν)-th moment of relaxationtime. For diffusion data, the kernel in eqn. (1) is e^(−Dt). Using eqn.(10) with ω=1, the negative first moment of diffusion and viscosity areobtained.

In addition to directly providing the parameters of interest, thesepulse sequences are optimized to estimate any desired moments ofrelaxation time and/or diffusion. In comparison to pulse sequences wherethe data is uniformly acquired in time, these pulse sequences can beused to acquire the data faster. For example, given ω>0, the error-barin the ω-th moment of T₁ relaxation is,

$\begin{matrix}{\left. \sigma_{\langle T_{1}^{\omega}\rangle} \right|_{{Uniform}\mspace{14mu}{Sampling}} = {\frac{\sigma_{\epsilon}}{{\Gamma(\omega)}\phi}t_{E}^{\omega}\sqrt{\sum\limits_{n = 1}^{N_{e}}n^{{2\;\omega} - 2}}}} & (23)\end{matrix}$Here N_(e) denotes the total number of echoes obtained from uniformsampling between a given lower and upper limit of relaxation time andt_(E) corresponds to the spacing between the echoes. In the case ofnon-uniform sampling, the error-bar in the ω-th moment is,

$\begin{matrix}{\left. \sigma_{\langle T_{1}^{\omega}\rangle} \right|_{{Non}\text{-}{Uniform}\mspace{14mu}{Sampling}} = {\frac{\sigma_{\epsilon}}{{\Gamma\left( {\omega + 1} \right)}\phi}\tau_{E}\sqrt{\; N_{e}}}} & (24)\end{matrix}$Here τ_(E) corresponds to the spacing between echoes in the τ domain.FIG. 5 is a plot showing the ratio of the error-bar in the momentsobtained from uniform sampling to non-uniform sampling. There are twopoints to note from this curve 510. First, there is a significantadvantage to non-uniform sampling as the error-bar in the computation ofmoments is lower (or at worst equal when ω=1). Second, from eqns. (23)and (24), the standard deviation of noise for uniform sampling has to besmaller than that for non-uniform sampling to obtain a desired error-barin the moments. Since the magnitude of the noise standard deviation inthe data is inversely proportional to the square-root of the acquisitiontime, to achieve the same uncertainty in the parameter the data fornon-uniform sampling can be acquired faster than that for uniformsampling.

Thus, according to some embodiments, a modified pulse sequence isprovided that can directly provide moments of relaxation-time ordiffusion distributions. This pulse sequence can be adapted to thedesired moment of relaxation-time or diffusion coefficient. The datafrom this pulse sequence provides direct estimates of fluid propertiessuch as average chain length and viscosity of a hydrocarbon. In general,in comparison with data acquired uniformly in time, this pulse sequenceis faster and has a lower error bar in computing the fluid properties.

A general expression for the sampling rate of non-uniformly spaced pulsesequence will now be described in further detail, according to someembodiments. When ω is not a negative integer,

T₁ ^(ω)

can be obtained by sampling uniformly in the τ domain, whereτ=t ^(ω)  (25)

Proof: For ω>0 and 1<ω<0, the proof has already been demonstrated ineqn. (11) and eqns. (12)-(17) respectively. When −2<ω≤1, from eqn. (8),n=2 and μ=ω+2. In this case, from eqn. (7),

$\begin{matrix}\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{1}{\phi\;{\Gamma(\mu)}}{\int_{0}^{\infty}{t^{\mu - 1}\frac{d^{2}}{{dt}^{2}}{M(t)}{dt}}}}} \\{= {\frac{1}{\phi\;{\Gamma(\mu)}}{\int_{0}^{\infty}{t^{\mu - 1}{\frac{d^{2}}{{dt}^{2}}\left\lbrack {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\rbrack}{dt}}}}} \\{= {{\frac{1}{\phi\;{\Gamma(\mu)}}\left\lbrack {t^{\mu - 1}\frac{d}{dt}\left\{ {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\}} \right\rbrack}_{t = 0}^{t = \infty} -}} \\{\frac{\left( {\mu - 1} \right)}{{\phi\Gamma}(\mu)}{\int_{0}^{\infty}{t^{\mu - 2}{\frac{d}{dt}\left\lbrack {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\rbrack}{dt}}}} \\{= {{{- \frac{1}{\phi}}{M^{(1)}(0)}\delta_{\mu\; 1}} - \frac{1}{{\phi\Gamma}\left( {\mu - 1} \right)}}} \\{\left\lbrack {t^{\mu - 2}\left\{ {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\}} \right\rbrack_{t = 0}^{t = \infty} +} \\{\frac{\left( {\mu - 2} \right)}{{\phi\Gamma}\left( {\mu - 1} \right)}{\int_{0}^{\infty}{{t^{\mu - 3}\left\lbrack {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\rbrack}{dt}}}} \\{= {{{- \frac{1}{\phi}}{M^{(1)}(0)}\delta_{\mu\; 1}} + {\frac{\left( {\mu - 2} \right)}{{\phi\Gamma}\left( {\mu - 1} \right)}{\int_{0}^{\infty}t^{\mu - 3}}}}} \\{\left\lbrack {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\rbrack{dt}} \\{= {{{- \frac{1}{\phi}}{M^{(1)}(0)}\delta_{\mu\; 1}} + {\frac{1}{{\phi\Gamma}(\omega)}{\int_{0}^{\infty}t^{\omega - 1}}}}} \\{\left\lbrack {{M(t)} - {M(0)} - {{tM}^{(1)}(0)}} \right\rbrack{dt}}\end{matrix} & (26)\end{matrix}$where δ_(μ1)=1 only when μ=1 and is zero otherwise and

${M^{1}(t)} = {\frac{{dM}(t)}{dt}.}$When μ=1 (corresponding to ω=1), the first term in eqn. (26) correspondsto the first derivative of the data at t=0 and is consistent with thedefinition of the negative first moment. In this case, the second termis zero. When −2<μ<1, the first term in eqn. (26) is zero. The secondterm is integrable at t=0. Using eqn. (25), the area under the integrand[M(t)−M(0)−tM⁽¹⁾(0)] can be used to directly estimate the ω-th moment.

In general, when n≥1, eqn. (7) can be written as

$\begin{matrix}\begin{matrix}{\left\langle T_{1}^{\omega} \right\rangle = {\frac{\left( {- 1} \right)}{{\phi\Gamma}(\mu)}{\int_{0}^{\infty}{t^{\mu - 1}\frac{d^{n}}{{dt}^{n}}{M(t)}{dt}}}}} \\{= {\frac{\left( {- 1} \right)^{n}}{{\phi\Gamma}(\mu)}{\int_{0}^{\infty}{t^{\mu - 1}{M^{(n)}(t)}{dt}}}}} \\{= {{\frac{\left( {- 1} \right)^{n - 1}}{\phi}{M^{({n - 1})}(0)}\delta_{\mu\; 1}} + {\frac{1}{{\phi\Gamma}\left( {\mu - n} \right)}{\int_{0}^{\infty}t^{\mu - n - 1}}}}} \\{\left\lbrack {{M(t)} - {\sum\limits_{m = 0}^{n - 1}{\frac{t^{m}}{m!}{M^{(m)}(0)}}}} \right\rbrack{dt}} \\{= {{\frac{\left( {- 1} \right)^{n - 1}}{\phi}{M^{({n - 1})}(0)}\delta_{\mu\; 1}} + {\frac{1}{{\phi\Gamma}(\omega)}{\int_{0}^{\infty}t^{\omega - 1}}}}} \\{\left\lbrack {{M(t)} - {\sum\limits_{m = 0}^{n - 1}{\frac{t^{m}}{m!}{M^{(m)}(0)}}}} \right\rbrack{dt}}\end{matrix} & (27)\end{matrix}$where

$\begin{matrix}{{M^{(m)}(t)} = {\frac{d^{m}}{{dt}^{m}}{M(t)}}} & (28)\end{matrix}$andM ⁽⁰⁾(t)=M(t)  (29)Note that:

$\begin{matrix}\left. {{M(t)} - {\sum\limits_{m = 0}^{n - 1}{\frac{t^{m}}{m!}{M^{(m)}(0)}}}}\rightarrow\left. {t^{n}\mspace{14mu}{as}\mspace{14mu} t}\rightarrow 0 \right. \right. & (30)\end{matrix}$

As before, the first term in eqn. (27) is valid only for negativeinteger values of ω and corresponds to the |ω|-th derivative of thedata. For negative values of ω where n<ω<(n−1), the first term in eqn.(27) is zero. Using eqn. (25), the area under the integrand

$\left\lbrack {{M(t)} - {\sum\limits_{m = 0}^{n - 1}{\frac{t^{m}}{m!}{M^{(m)}(0)}}}} \right\rbrack$in the τ domain directly provides the ω-th moment of relaxation time.

FIG. 6A shows an NMR tool being deployed in a wellbore, according toembodiments. Wireline truck 610 is deploying wireline cable 612 intowell 630 via well head 620. Wireline tool 640 is disposed on the end ofthe cable 612 in a subterranean formation 600. Tool 640 includes an NMRtool, such as Schlumberger's Combinable Magnetic Resonance (CMR) Tool,or Schlumberger's MR Scanner Tool. According to some embodiments the NMRtool is combined with other downhole tools such as Schlumberger'sModular Formation Dynamics Tester downhole sampling tool. According tosome embodiments, the NMR tool uses a pulse sequence as describedherein. Data from the NMR tool can be recorded and/or processed downholewithin tool 640, and/or can be transmitted to wireline truck 610 forrecording and/or processing. According to embodiments, the NMR tool canbe controlled, including the selection of pulse sequences, locally withprocessing systems within tool 640 and/or from the surface usingprocessing systems for example in wireline truck 610.

FIG. 6B is a schematic of a processing system used to control, recordand process data from an NMR tool. Processing system 650 includes one ormore central processing units 674, storage system 672, communicationsand input/output modules 670, a user display 676 and a user input system678. Input/output modules 670 include modules to communicate with andcontrol the NMR tool such that modified pulse sequences can be used asdescribed herein.

FIG. 7 illustrates another example of a wellsite system in which an NMRtool can be employed, according to some embodiments. The wellsite can beonshore or offshore. In this exemplary system, a borehole 711 is formedin subsurface formations by rotary drilling in a manner that is wellknown. Embodiments of the invention can also use directional drilling,as will be described hereinafter.

A drill string 712 is suspended within the borehole 711 and has a bottomhole assembly 700 which includes a drill bit 705 at its lower end. Thesurface system includes platform and derrick assembly 710 positionedover the borehole 711, the assembly 710 including a rotary table 716,kelly 717, hook 718 and rotary swivel 719. The drill string 712 isrotated by the rotary table 716, energized by means not shown, whichengages the kelly 717 at the upper end of the drill string. The drillstring 712 is suspended from a hook 718, attached to a traveling block(also not shown), through the kelly 717 and a rotary swivel 719 whichpermits rotation of the drill string relative to the hook. As is wellknown, a top drive system could alternatively be used.

In the example of this embodiment, the surface system further includesdrilling fluid or mud 726 stored in a pit 727 formed at the well site. Apump 729 delivers the drilling fluid 726 to the interior of the drillstring 712 via a port in the swivel 719, causing the drilling fluid toflow downwardly through the drill string 712 as indicated by thedirectional arrow 708. The drilling fluid exits the drill string 712 viaports in the drill bit 705, and then circulates upwardly through theannulus region between the outside of the drill string and the wall ofthe borehole, as indicated by the directional arrows 709. In this wellknown manner, the drilling fluid lubricates the drill bit 705 andcarries formation cuttings up to the surface as it is returned to thepit 727 for recirculation.

The bottom hole assembly 700 of the illustrated embodiment alogging-while-drilling (LWD) module 720, a measuring-while-drilling(MWD) module 730, a roto-steerable system and motor, and drill bit 705.

The LWD module 720 is housed in a special type of drill collar, as isknown in the art, and can contain one or a plurality of known types oflogging tools. It will also be understood that more than one LWD and/orMWD module can be employed, e.g. as represented at 720A. (References,throughout, to a module at the position of 720 can alternatively mean amodule at the position of 720A as well.) The LWD module includescapabilities for measuring, processing, and storing information, as wellas for communicating with the surface equipment. In the presentembodiment, the LWD module includes an NMR tool such as Schlumberger'sProVISION NMR tool.

The MWD module 730 is also housed in a special type of drill collar, asis known in the art, and can contain one or more devices for measuringcharacteristics of the drill string and drill bit. The MWD tool furtherincludes an apparatus (not shown) for generating electrical power to thedownhole system. This may typically include a mud turbine generatorpowered by the flow of the drilling fluid, it being understood thatother power and/or battery systems may be employed. In the presentembodiment, the MWD module includes one or more of the following typesof measuring devices: a weight-on-bit measuring device, a torquemeasuring device, a vibration measuring device, a shock measuringdevice, a stick slip measuring device, a direction measuring device, andan inclination measuring device.

According to some embodiments, the techniques described herein areapplied to estimate properties from other types of excitations than NMR.For example the techniques can be applied to other types of data whereexcitations induce an exponential decay. According to some embodiments,the techniques are applied to inducing a non-uniform sequence ofpressure changes in a material. According to some other embodiments, thetechniques are applied to induce a non-uniform sequence of temperaturechanges in a material.

While many examples have been described herein with respect to designinginterval spacings in the time domain, according to some embodiments, thetechniques described herein are applied to interval spacings in otherdomains, such as designing interval spacings in the frequency domain.

Whereas many alterations and modifications of the present disclosurewill no doubt become apparent to a person of ordinary skill in the artafter having read the foregoing description, it is to be understood thatthe particular embodiments shown and described by way of illustrationare in no way intended to be considered limiting. Further, thedisclosure has been described with reference to particular preferredembodiments, but variations within the spirit and scope of thedisclosure will occur to those skilled in the art. It is noted that theforegoing examples have been provided merely for the purpose ofexplanation and are in no way to be construed as limiting of the presentdisclosure. While the present disclosure has been described withreference to exemplary embodiments, it is understood that the words,which have been used herein, are words of description and illustration,rather than words of limitation. Changes may be made, within the purviewof the appended claims, as presently stated and as amended, withoutdeparting from the scope and spirit of the present disclosure in itsaspects. Although the present disclosure has been described herein withreference to particular means, materials and embodiments, the presentdisclosure is not intended to be limited to the particulars disclosedherein; rather, the present disclosure extends to all functionallyequivalent structures, methods and uses, such as are within the scope ofthe appended claims.

What is claimed is:
 1. A method of estimating a property of a materialusing nuclear magnetic resonance (NMR), the method comprising: inducinga measurable change in the material using an NMR pulse; after apredetermined time interval, making a measurement of a material responseto the induction using NMR, wherein the material response exponentiallydecays following the induction; repeating the inducing and the measuringat non-uniform predetermined time intervals, wherein the non-uniformpredetermined time intervals are determined based on the property to beestimated; and estimating the property based at least in part on themeasurements of the exponentially decaying material response.
 2. Amethod according to claim 1 wherein the NMR pulse is a transversemagnetic pulse.
 3. A method according to claim 1 wherein thepredetermined time intervals are non-uniform in both linear andlogarithmic scales.
 4. A method according to claim 1 wherein thenon-uniform predetermined time intervals follow a power law.
 5. A methodaccording to claim 1 wherein the non-uniform predetermined timeintervals allow for the property to be directly computed from themeasurements of the exponentially decaying material response.
 6. Amethod according to claim 1 wherein the measurements are performed aspart of one or more oilfield services.
 7. A method according to claim 6wherein the measurements are made downhole.
 8. A method according toclaim 1 wherein the material is a fluid.
 9. A method according to claim8 wherein the property is a fluid property.
 10. A method according toclaim 1 wherein the material is one or more rock solids and the propertyis a petrophysical property of the one or more rock solids.
 11. A methodaccording to claim 1 wherein the non-uniform predetermined timeintervals are determined based on a moment of relaxation-time ordiffusion distribution.
 12. A method according to claim 1 wherein thematerial is a mixture of solid and liquid.
 13. A method according toclaim 12 wherein the property comprises a plurality of properties andthe plurality of properties comprises a fluid property and apetrophysical property in porous solids.
 14. A nuclear magneticresonance (NMR) system for estimating a property of a material, the NMRsystem comprising: a NMR tool adapted to (i) induce a measurable changein the material using an NMR pulse and (ii) measure a material responseto the induction after a predetermined time interval; and a processingsystem adapted to (i) repeat the induction and the measurement atnon-uniform predetermined time intervals, wherein the non-uniformpredetermined time intervals are determined based on the property to beestimated, and (ii) estimate the property based at least in part on themeasurements.
 15. A system according to claim 14 wherein the NMR tooltransmits a transverse magnetic pulse using an antenna.
 16. A systemaccording to claim 14 wherein the non-uniform predetermined timeintervals are non-uniform in both linear and logarithmic scales.
 17. Asystem according to claim 14 wherein the non-uniform predetermined timeintervals follow a power law.
 18. A system according to claim 14 whereinthe non-uniform predetermined time intervals allow for the property tobe directly computed from the measurements.
 19. A system according toclaim 14 wherein the NMR system is part of an oilfield services tool.20. A system according to claim 14 further comprising a downholedeployable tool body that houses the NMR tool.
 21. A system according toclaim 14 wherein the material is a fluid.
 22. A system according toclaim 21 wherein the property is a fluid property.
 23. A systemaccording to claim 14 wherein the material is one or more rock solidsand the property is a petrophysical property of the one or more rocksolids.
 24. A system according to claim 14 wherein the non-uniformpredetermined time intervals are determined based on a moment ofrelaxation-time or diffusion distribution.
 25. A system according toclaim 14 wherein the material is a mixture of solid and liquid.
 26. Amethod according to claim 25 wherein the property is a fluid property inporous solids.